{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def function_f(x):\n",
    "    return 1.0\n",
    "def function_u(x):\n",
    "    return -x*x/2\n",
    "\n",
    "def basic_function(i,x):\n",
    "    if i== 0:\n",
    "        return (x-1)*x*1/2\n",
    "    elif i==1:\n",
    "        return (x+1)*(x-1)*-1\n",
    "    elif i==2:\n",
    "        return (x+1)*x*1/2\n",
    "\n",
    "def basic_function_grad(i,x):\n",
    "    if i ==0:\n",
    "        return x-0.5\n",
    "    elif i==1:\n",
    "        return -2*x\n",
    "    elif i==2:\n",
    "        return x+0.5\n",
    "\n",
    "def dealwithboundary(A,F,x_l,x_r):\n",
    "    u0 = function_u(x_l)\n",
    "    u1 = function_u(x_r)\n",
    "    size = len(F)\n",
    "    for i in range(1,size-1):\n",
    "        F[i] -= A[i][0] * u0 + A[i][size-1] *u1\n",
    "        A[i][0] = 0\n",
    "        A[i][size -1] = 0\n",
    "        A[0][i] = 0\n",
    "        A[size-1][i] = 0 \n",
    "    F[0] = A[0][0] * u0\n",
    "    F[size-1] = A[size-1][size-1]*u1\n",
    "    return \n",
    "\n",
    "Gauss_Pnt_p = [-0.9061798459386640, -0.5384693101056830, 0.0000000000000000, 0.5384693101056830, 0.9061798459386640]\n",
    "Gauss_Pnt_w = [0.2369268850561890, 0.4786286704993660, 0.5688888888888880, 0.4786286704993660, 0.2369268850561890]\n",
    "\n",
    "x_l = 0.0\n",
    "x_r = 1.0\n",
    "n_size = 10\n",
    "h = (x_r - x_l) / n_size\n",
    "\n",
    "\n",
    "def fem_2grad(n_size):\n",
    "    A = np.zeros((2*n_size+1, 2*n_size+1))\n",
    "    F = np.zeros((2*n_size+1))\n",
    "    for i in range(n_size):\n",
    "        x_lhs = i*h + x_l\n",
    "        #########拼装刚度矩阵\n",
    "        for j in range(3):\n",
    "            for k in range(3):\n",
    "                for l in range(5):\n",
    "                    A[2*i+j][2*i+k] += basic_function_grad(j,Gauss_Pnt_p[l])*basic_function_grad(k,Gauss_Pnt_p[l])*Gauss_Pnt_w[l]*2/h\n",
    "            \n",
    "            for l in range(5):\n",
    "                F[2*i+j] += Gauss_Pnt_w[l] * basic_function(j, Gauss_Pnt_p[l]) * function_f((Gauss_Pnt_p[l]+1)/2 *h + x_lhs) * h /2\n",
    "    return A,F\n",
    "\n",
    "A_fine, F_fine = fem_2grad(2*n_size)\n",
    "\n",
    "A_coarse, F_coarse = fem_2grad(n_size)\n",
    "\n",
    "dealwithboundary(A_fine,F_fine,x_l,x_r)\n",
    "dealwithboundary(A_coarse,F_coarse,x_l,x_r)\n",
    "\n",
    "u_fine = np.linalg.solve(A_fine,F_fine)\n",
    "u_coarse = np.linalg.solve(A_coarse,F_coarse)"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "[[ 23.33333333 -26.66666667   3.33333333   0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [-26.66666667  53.33333333 -26.66666667   0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  3.33333333 -26.66666667  46.66666667 -26.66666667   3.33333333\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.         -26.66666667  53.33333333 -26.66666667\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           3.33333333 -26.66666667  46.66666667\n",
      "  -26.66666667   3.33333333   0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.         -26.66666667\n",
      "   53.33333333 -26.66666667   0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           3.33333333\n",
      "  -26.66666667  46.66666667 -26.66666667   3.33333333   0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.         -26.66666667  53.33333333 -26.66666667   0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           3.33333333 -26.66666667  46.66666667 -26.66666667\n",
      "    3.33333333   0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.         -26.66666667  53.33333333\n",
      "  -26.66666667   0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           3.33333333 -26.66666667\n",
      "   46.66666667 -26.66666667   3.33333333   0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "  -26.66666667  53.33333333 -26.66666667   0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    3.33333333 -26.66666667  46.66666667 -26.66666667   3.33333333\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.         -26.66666667  53.33333333 -26.66666667\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           3.33333333 -26.66666667  46.66666667\n",
      "  -26.66666667   3.33333333   0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.         -26.66666667\n",
      "   53.33333333 -26.66666667   0.           0.           0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           3.33333333\n",
      "  -26.66666667  46.66666667 -26.66666667   3.33333333   0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.         -26.66666667  53.33333333 -26.66666667   0.\n",
      "    0.        ]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           3.33333333 -26.66666667  46.66666667 -26.66666667\n",
      "    3.33333333]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.         -26.66666667  53.33333333\n",
      "  -26.66666667]\n",
      " [  0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           0.           0.\n",
      "    0.           0.           0.           3.33333333 -26.66666667\n",
      "   23.33333333]]\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "source": [
    "print(u_fine)\n"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "[ 0.      -0.00125 -0.005   -0.01125 -0.02    -0.03125 -0.045   -0.06125\n",
      " -0.08    -0.10125 -0.125   -0.15125 -0.18    -0.21125 -0.245   -0.28125\n",
      " -0.32    -0.36125 -0.405   -0.45125 -0.5    ]\n"
     ]
    }
   ],
   "metadata": {}
  }
 ],
 "metadata": {
  "orig_nbformat": 4,
  "language_info": {
   "name": "python",
   "version": "3.7.10",
   "mimetype": "text/x-python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "pygments_lexer": "ipython3",
   "nbconvert_exporter": "python",
   "file_extension": ".py"
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.7.10 64-bit ('base': conda)"
  },
  "interpreter": {
   "hash": "2c21737055494eda6dd7820cee4b5069ee6eef9e680dd0a23b71fc1db8cdc626"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}